# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# GIGER / LANZ / DEVRIES 2018

# PSRM-OA-2017-0118 -- REPLICATION APPENDIX

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

rm(list=ls(all=TRUE)) 

library(foreign)
library(MASS)  
library(sandwich)	
library(zoo)
library(lmtest)
library(Formula)
library(betareg)
library(mfx)
library(lme4)
library(mlmRev)
library(arm)
library(stargazer)
library(RColorBrewer)
library(nloptr)
library(stringr)
library(readstata13)
library(reshape2)
library(ggplot2)
library(randomizr)
library(randomizr)

load("DATA_PSRM-OA-2017-0118.RData")

data_comp      <- data                           # complete dataset with all candidates (contacted and not contacted)
data           <- subset(data, answer!= "NA")    # data with contacted candidates (main data for analysis)
data_awn       <- subset(data, answer== "yes")   # data with contacted candidates & anwsers
data_ninc      <- subset(data, inc_r==0)         # data without incumbents
data_valais    <- subset(data, sender!=23)       # data without canton of Valais
data_staff     <- subset(data, staff==0)         # data without candidates with payed staff
data_list      <- subset(data, list_alph==0)     # data wehere listpos is non-alphabetical
data_list_alph <- subset(data, list_alph==1)     # data wehere listpos is alphabetical

options(warn=-1)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A1 Balance Tests

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A1: Candidate characteristics by treatment
# ------------------------------------------------------------------------------------------- #

# party treatment (first variable from top)

table(data$treatment_party)
table(data$treatment_cant, data$treatment_party)
prop.table(table(data$treatment_cant, data$treatment_party), 1)

t.test(data$treatment_party_r[data$treatment_cant == "out canton"],
       data$treatment_party_r[data$treatment_cant == "in canton"])

# safety (second variable from top)

mean(data$safety[data$treatment_cant == "out canton"])
sd(data$safety[data$treatment_cant == "out canton"])
mean(data$safety[data$treatment_cant == "in canton"])
sd(data$safety[data$treatment_cant == "in canton"])

t.test(data$safety[data$treatment_cant == "out canton"],
       data$safety[data$treatment_cant == "in canton"])

# age (third variable from top)

mean(data$age[data$treatment_cant == "out canton"])
sd(data$age[data$treatment_cant == "out canton"])
mean(data$age[data$treatment_cant == "in canton"])
sd(data$age[data$treatment_cant == "in canton"])

t.test(data$age[data$treatment_cant == "out canton"],
       data$age[data$treatment_cant == "in canton"])

# gender (fourth variable from top)

table(data$male)
table(data$treatment_cant, data$male)

prop.table(table(data$treatment_cant, data$male), 1)

t.test(data$male_r[data$treatment_cant == "out canton"],
       data$male_r[data$treatment_cant == "in canton"])

# party (fifth variable from top)

table(data$party1)
table(data$treatment_cant, data$party1)
prop.table(table(data$treatment_cant, data$party1), 1)

t.test(data$party1[data$treatment_cant == "out canton"],
       data$party1[data$treatment_cant == "in canton"])

table(data$party2)
table(data$treatment_cant, data$party2)
prop.table(table(data$treatment_cant, data$party2), 1)

t.test(data$party2[data$treatment_cant == "out canton"],
       data$party2[data$treatment_cant == "in canton"])

table(data$party3)
table(data$treatment_cant, data$party3)
prop.table(table(data$treatment_cant, data$party3), 1)

t.test(data$party3[data$treatment_cant == "out canton"],
       data$party3[data$treatment_cant == "in canton"])

table(data$party4)
table(data$treatment_cant, data$party4)
prop.table(table(data$treatment_cant, data$party4), 1)

t.test(data$party4[data$treatment_cant == "out canton"],
       data$party4[data$treatment_cant == "in canton"])

table(data$party5)
table(data$treatment_cant, data$party5)
prop.table(table(data$treatment_cant, data$party5), 1)

t.test(data$party5[data$treatment_cant == "out canton"],
       data$party5[data$treatment_cant == "in canton"])

# language (sixth variable from top)

table(data$lang)
table(data$treatment_cant, data$lang)
prop.table(table(data$treatment_cant, data$lang), 1)

t.test(data$lang_r[data$treatment_cant == "out canton"],
       data$lang_r[data$treatment_cant == "in canton"])

# nr seats (seventh variable from top)

mean(data$seats[data$treatment_cant == "out canton"])
sd(data$seats[data$treatment_cant == "out canton"])
mean(data$seats[data$treatment_cant == "in canton"])
sd(data$seats[data$treatment_cant == "in canton"])

t.test(data$seats[data$treatment_cant == "out canton"],
       data$seats[data$treatment_cant == "in canton"])


# Table A2: RI logistic regression (outcome: treatment in-canton)
# ------------------------------------------------------------------------------------------- #

# m0

M0app <- glmer(treatment_cant ~
               + (1 | canton), data=data, 
               family = binomial(link = "logit"),  
               control=glmerControl(optimizer="nloptwrap"))

summary(M0app)   

# m1

M1app <- glmer(treatment_cant ~
               + treatment_party_r
               + (1 | canton), data=data, 
               family = binomial(link = "logit"),  
               control=glmerControl(optimizer="nloptwrap"))

summary(M1app)   

# m2

M2app <- glmer(treatment_cant ~ 
               + treatment_party_r
               + age 
               + male_r
               + party1 + party2 + party4 + party5 
               + lang_r
               + (1 | canton), data=data, 
               family = binomial(link = "logit"),  
               control=glmerControl(optimizer="nloptwrap"))

summary(M2app)    

# m3

M3app <- glmer(treatment_cant ~
               + treatment_party_r
               + safety
               + age 
               + male_r
               + party1 + party2 + party4 + party5 
               + lang_r
               + (1 | canton), data=data, 
               family = binomial(link = "logit"),  
               control=glmerControl(optimizer="nloptwrap"))

summary(M3app)    

# Table A3: Likelihood ratio test (M1, M2, M3 vs. M0)
# ------------------------------------------------------------------------------------------- #

lrtest(M0app, M1app)
lrtest(M0app, M2app)
lrtest(M0app, M3app)

# clean up

rm(M0app, M1app, M2app, M3app)

# Figure A1: Density plot of simulated log likelihoods and the log likelihood of the actual 
# treatment (M3 in Table A6, purple line) (export: 12in X 4.5in)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,1))
par(mar=c(3,1.5,3,1.5), oma=c(0,0,0,0))

ll <- data.frame()

for(i in 1:5000){
  treat <- complete_ra(N=nrow(data))
  data["treat"] <- treat
  
  Msim <- glmer(treat ~ 
                     + treatment_party_r
                     + safety
                     + age 
                     + male_r
                     + party1 + party2 + party4 + party5 
                     + lang_r
                     + (1 | canton), data=data, 
                     family = binomial(link = "logit"),  
                     control=glmerControl(optimizer="nloptwrap"))
  
  ll <- rbind(ll, c(logLik(Msim)))
}

colnames(ll)<- c("loglike")

Mtreat <- glmer(treatment_cant_r ~ 
                 + treatment_party_r
                 + safety
                 + age 
                 + male_r
                 + party1 + party2 + party4 + party5 
                 + lang_r
                 + (1 | canton), data=data, 
                 family = binomial(link = "logit"),  
                 control=glmerControl(optimizer="nloptwrap"))

loglik <- logLik(Mtreat)

test <- ll$loglike > loglik
table(test)                                   
pvalue <- length(test[test==TRUE])/length(test[test==FALSE]) 
pvalue

dens_ll <- density(ll$loglike)

plot(dens_ll, ylab=" ", yaxt="n", xlim=c(-460,-440), ylim=c(0,0.26),  
     main="Balance test treatment", type="n", 
     xaxs='i',yaxs='i', cex.main=0.85, cex.axis=0.85)

abline(v=seq(-455,-435,by=5), col="grey80")
abline(h=seq(0,0.27,by=0.05), col="grey80")
abline(h=0, col="black")
points(dens_ll, type="l", lwd= 2)

mean(ll$loglike)
loglik
text(-440, 0.24, "N simulations = 5,000; Mean ll = -453", cex=0.85, pos=2)

abline(v=seq(0,-450 , by=-450), col="purple4", lwd=2)
text(-450, 0.16, "ll actual treatment = -450", cex=0.85, pos=4)

# save chart

dev.copy(pdf,"figureA1.pdf", width=12, height=4.5)
dev.off()

# clean up

rm(dens_ll, loglik, Msim, Mtreat, pvalue, test, treat, ll, i)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A2 Composition of the Sample

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A4: Candidate characteristics by treatment
# ------------------------------------------------------------------------------------------- #

# age (first variable from top)

mean(data_comp$age[data_comp$contacted == "no"], na.rm=TRUE)
sd(data_comp$age[data_comp$contacted == "no"], na.rm=TRUE)

mean(data_comp$age[data_comp$contacted == "yes"], na.rm=TRUE)
sd(data_comp$age[data_comp$contacted == "yes"], na.rm=TRUE)

t.test(data_comp$age[data_comp$contacted == "no"],
       data_comp$age[data_comp$contacted == "yes"])

# gender (second varible from top)

table(data_comp$male)
table(data_comp$contacted, data_comp$male)
prop.table(table(data_comp$contacted, data_comp$male), 1)

t.test(data_comp$male_r[data_comp$contacted == "no"],
       data_comp$male_r[data_comp$contacted == "yes"])

# party (third variable from top)

table(data_comp$party1)
table(data_comp$contacted, data_comp$party1)
prop.table(table(data_comp$contacted, data_comp$party1), 1)

t.test(data_comp$party1[data_comp$contacted == "no"],
       data_comp$party1[data_comp$contacted == "yes"])

table(data_comp$party2)
table(data_comp$contacted, data_comp$party2)
prop.table(table(data_comp$contacted, data_comp$party2), 1)

t.test(data_comp$party2[data_comp$contacted == "no"],
       data_comp$party2[data_comp$contacted == "yes"])

table(data_comp$party3)
table(data_comp$contacted, data_comp$party3)
prop.table(table(data_comp$contacted, data_comp$party3), 1)

t.test(data_comp$party3[data_comp$contacted == "no"],
       data_comp$party3[data_comp$contacted == "yes"])

table(data_comp$party4)
table(data_comp$contacted, data_comp$party4)
prop.table(table(data_comp$contacted, data_comp$party4), 1)

t.test(data_comp$party4[data_comp$contacted == "no"],
       data_comp$party4[data_comp$contacted == "yes"])

table(data_comp$party5)
table(data_comp$contacted, data_comp$party5)
prop.table(table(data_comp$contacted, data_comp$party5), 1)

t.test(data_comp$party5[data_comp$contacted == "no"],
       data_comp$party5[data_comp$contacted == "yes"])

# language (fourth variable from top)

table(data_comp$lang)
table(data_comp$contacted, data_comp$lang)
prop.table(table(data_comp$contacted, data_comp$lang), 1)

t.test(data_comp$lang_r[data_comp$contacted == "no"],
       data_comp$lang_r[data_comp$contacted == "yes"])

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A3 Descriptive Statistics: Electoral safety

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Figure A2: Density plot: electoral safety (export: 12in X 4.5in)
# ------------------------------------------------------------------------------------------- #

dev.off()
par(mfrow=c(1,1))
par(mar=c(3,1.5,3,1.5), oma=c(0,0,0,0))

dens_saf <- density(data$safety)

plot(dens_saf, ylab=" ", yaxt="n",  xlim=c(0,200), ylim=c(0,0.02), 
     main=NA, type="n" , xaxs='i',yaxs='i', cex.main=0.85,  cex.axis=0.85)

abline(v=seq(0,175,by=50), col="grey80")
abline(h=seq(0,0.015,by=0.005), col="grey80")
abline(h=0, col="black")
points(dens_saf, type="l", lwd= 2)

mean(data$safety)
length(data$safety)

text(160, 0.018, "N=660; Mean=62.4", cex=0.9)

# save chart

dev.copy(pdf,"figureA2.pdf", width=12, height=4.5)
dev.off()

# clean up

rm(dens_saf)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A4 Descriptive Statistics: Reported Importance of Constituency Service

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A5: Importance of constituency service and responsiveness
# ------------------------------------------------------------------------------------------- 

# first/second column (answer no)

table(data$service, data$answer)
prop.table(table(data$service, data$answer), 1)

# third column (all)

table(data$service)
prop.table(table(data$service))

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A5 Treatment Effect: Further Analyses

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A6: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

# m1

M1 <- glmer(answer ~ 
            + treatment_cant_r
            + treatment_party_r
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1)

# m2

M2 <- glmer(answer ~ 
            + treatment_cant_r
            + treatment_party_r
            + age
            + male_r
            + party1 + party2 + party4 + party5 
            + lang_r
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M2)

# m3

M3 <- glmer(answer ~ 
            + treatment_cant_r
            + treatment_party_r
            + safety
            + age
            + male_r
            + party1 + party2 + party4 + party5 
            + lang_r
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M3)

# clean up

rm(M1, M2, M3)


# Figure A3: Responsiveness split by cantonal treatment (M1 is also in rep_psrm_main.R)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,1))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

M1 <- glmer(answer ~ 
              treatment_cant_r
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1)

m.beta <- fixef(M1)
v.beta <- vcov(M1)

XX0a <- cbind(1, 
              0)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.star0a <- XX0a %*% t(BETA)
p.hat0a <- 1/(1+exp(-y.star0a))

XX0b <- cbind(1, 
              1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.star0b <- XX0b %*% t(BETA)
p.hat0b <- 1/(1+exp(-y.star0b))

mean(p.hat0a) # (answer prob if treatment = OUT canton)
mean(p.hat0b) # (answer prob if treatment = IN canton)

fd1 <- p.hat0b - p.hat0a
sort(fd1)[c(37.5,1462.5)]

mean(fd1) # (fd in canton vs. vs canton) 

pred=matrix(0,2,3)

pred[1,1]=apply(p.hat0a,1,FUN=median) # out canton
pred[1,2]=apply(p.hat0a,1,FUN=quantile,probs = .025)
pred[1,3]=apply(p.hat0a,1,FUN=quantile,probs = 0.975)

pred[2,1]=apply(p.hat0b,1,FUN=median) # in canton
pred[2,2]=apply(p.hat0b,1,FUN=quantile,probs = .025)
pred[2,3]=apply(p.hat0b,1,FUN=quantile,probs = 0.975)

plot (data$answer_dich, data$answer_dich, 
      xlim=c(0.5,2.5), ylim=c(0.4,0.9),
      ylab="Pr (Answer)", xlab="Treatment", 
      type="n", xaxt="n", yaxs="i", 
      mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)

abline(h=seq(0.5,0.8,by=0.1), col="grey92")

segments(c(1:2), pred[,2], c(1:2), pred[,3], lwd= 2 , col="black")
points(c(1:2),pred[,1], pch=15, cex=0.7 , col="black")

text(1, 0.86, "Out canton", pos=1, cex=0.75, col="black")
text(2, 0.86, "In canton", pos=1, cex=0.75, col="black")

# save chart

dev.copy(pdf,"figureA3.pdf", width=5.5, height=4)
dev.off()

# clean up

rm(BETA, fd1, p.hat0a, p.hat0b, pred, XX0a, XX0b, y.star0a, y.star0b, m.beta, M1, v.beta)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A6 Visualizations M2, and M3, Table 2

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Figure A4: Treatments and responsiveness (M2 is also in rep_psrm_main.R; export 9in X 4 in)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

M2 <- glmer(answer ~ 
            + treatment_cant_r
            + safety
            + I(treatment_cant_r*safety)
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M2)

m.beta <- fixef(M2)
v.beta <- vcov(M2)

# left panel

XX0 <- cbind(1, 
             0, 
             seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), 
             0*seq(min(data$safety[!is.na(data$safety)]), 150, by = 1))

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             1, 
             seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), 
             1*seq(min(data$safety[!is.na(data$safety)]), 150, by = 1))

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety[!is.na(data$safety)]),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data$safety[!is.na(data$safety)]),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data$safety[data$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety[!is.na(data$safety)]),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data$safety[!is.na(data$safety)]),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA4.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, 
   p.hat0, p.hat1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M2)


# Figure A5: Treatments and responsiveness (M3 is also in rep_psrm_main.R; export 9in X 4 in)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

M3 <- glmer(answer ~ 
            + treatment_party_r
            + treatment_cant_r
            + safety
            + I(treatment_cant_r*safety)
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M3)

m.beta <- fixef(M3)
v.beta <- vcov(M3)

# left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), 
             0*seq(min(data$safety[!is.na(data$safety)]), 150, by = 1))

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), 
             1*seq(min(data$safety[!is.na(data$safety)]), 150, by = 1))

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

plot (seq(min(data$safety[!is.na(data$safety)]),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data$safety[!is.na(data$safety)]),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data$safety[data$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety[!is.na(data$safety)]),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data$safety[!is.na(data$safety)]),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA5.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, 
   p.hat0, p.hat1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M3)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A7 Robustness Check: List Position

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A7: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

data_list$list_pos <-(data_list$list_pos)/(data_list$seats)
data_list$list_pos <-(data_list$list_pos)-min(data_list$list_pos)
data_list$list_pos <-(data_list$list_pos)*(1/max(data_list$list_pos))
data_list$list_pos <-max(data_list$list_pos)+min(data_list$list_pos)-(data_list$list_pos)

M1pos <- glmer(answer ~ 
                  + treatment_party_r
                  + treatment_cant_r
                  + list_pos
                  + I(treatment_cant_r*list_pos)
                  + age 
                  + male_r
                  + party1 + party2 + party4 + party5 
                  + lang_r
                  + (1 | canton), data=data_list, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1pos)

# Figure A6: Responsiveness spit by cantonal treatment and across ballot position (export: 9 in X 4in)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

m.beta <- fixef(M1pos)
v.beta <- vcov(M1pos)

# left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data_list$list_pos), 1, by = 0.05), 
             0*seq(min(data_list$list_pos), 1, by = 0.05),
             mean(data_list$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data_list$list_pos), 1, by = 0.05), 
             1*seq(min(data_list$list_pos), 1, by = 0.05),
             mean(data_list$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data_list$list_pos),1, by = 0.05), p.ci.XX0[2,], xlab="Quality of ballot position", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data_list$list_pos),1),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,1, by = 0.2), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,1,by=0.2), col="grey92")

points(seq(min(data_list$list_pos), 1, by = 0.05), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data_list$list_pos), 1, by = 0.05), p.ci.XX1[2,], type="l", lwd= 2)

text(0.05, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(0.05, p.ci.XX1[2,1]-0.045, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data_list$list_pos, 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data_list$list_pos),1, by = 0.05), p.ci.XX0[2,], xlab="Quality of ballot position", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data_list$list_pos),1),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,1, by = 0.2), cex.axis=0.7)
abline(h=seq(-0.2,0.4,by=0.2), col="grey92")
abline(v=seq(0,1,by=0.2), col="grey92")
abline(h=0, col="black")

points(seq(min(data_list$list_pos), 1, by = 0.05), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data_list$list_pos), 1, by = 0.05), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data_list$list_pos), 1, by = 0.05), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA6.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, 
   p.hat0, p.hat1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M1pos)

# Table A8: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

data_list_alph$list_pos_alph <-(data_list_alph$list_pos)/(data_list_alph$seats)
data_list_alph$list_pos_alph <-(data_list_alph$list_pos_alph)-min(data_list_alph$list_pos_alph)
data_list_alph$list_pos_alph <-(data_list_alph$list_pos_alph)*(1/max(data_list_alph$list_pos_alph))
data_list_alph$list_pos_alph <-max(data_list_alph$list_pos_alph)+min(data_list_alph$list_pos_alph)-(data_list_alph$list_pos_alph)


M1alph <- glmer(answer ~ 
            + treatment_cant_r
            + list_pos_alph
            + I(treatment_cant_r*list_pos_alph)
            + (1 | canton), data=data_list_alph, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))

summary(M1alph)

# clean up

rm(M1alph)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A8 Robustness Check: Squared Electoral Safety

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A9: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

M1sq <- glmer(answer ~ 
                   + treatment_party_r
                   + treatment_cant_r
                   + safety
                   + I(treatment_cant_r*safety)
                   + I(safety*safety)
                   + I(treatment_cant_r*safety*safety)
                   + age 
                   + male_r
                   + party1 + party2 + party4 + party5 
                   + lang_r
                   + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1sq)

# Figure A7: Responsiveness spit by cantonal treatment and across squared electoral safety (export: 9 in X 4in)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

m.beta <- fixef(M1sq)
v.beta <- vcov(M1sq)

# left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data$safety), 150, by = 1), 
             0*seq(min(data$safety), 150, by = 1),
             seq(min(data$safety), 150, by = 1)*seq(min(data$safety), 150, by = 1),
             0*seq(min(data$safety), 150, by = 1)*seq(min(data$safety), 150, by = 1),
             mean(data$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data$safety), 150, by = 1), 
             1*seq(min(data$safety), 150, by = 1),
             seq(min(data$safety), 150, by = 1)*seq(min(data$safety), 150, by = 1),
             1*seq(min(data$safety), 150, by = 1)*seq(min(data$safety), 150, by = 1),
             mean(data$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data$safety),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data$safety), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data$safety), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data$safety[data$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data$safety),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data$safety), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data$safety), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA7.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, 
   p.hat0, p.hat1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M1sq)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A9 Robustness Check: Controlling for Incumbency Status

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A10: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

M1inc <- glmer(answer ~ 
                + treatment_party_r
                + treatment_cant_r
                + safety
                + I(treatment_cant_r*safety)
                + inc_r
                + age 
                + male_r
                + party1 + party2 + party4 + party5 
                + lang_r
                + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1inc)

# Figure A8: Responsiveness spit by cantonal treatment and across electoral safety (incumbency 
# included) (export: 9 in X 4in)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

m.beta <- fixef(M1inc)
v.beta <- vcov(M1inc)

# left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data$safety), 150, by = 1), 
             0*seq(min(data$safety), 150, by = 1),
             0,
             mean(data$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data$safety), 150, by = 1), 
             1*seq(min(data$safety), 150, by = 1),
             0,
             mean(data$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data$safety),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data$safety), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data$safety), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data$safety[data$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data$safety),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data$safety), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data$safety), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA8.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, 
   p.hat0, p.hat1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M1inc)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A10 Robustness Check: Incumbents Excluded from Analysis

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A11: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

M1ninc <- glmer(answer ~ 
                + treatment_party_r
                + treatment_cant_r
                + safety
                + I(treatment_cant_r*safety)
                + age 
                + male_r
                + party1 + party2 + party4 + party5 
                + lang_r
                + (1 | canton), data=data_ninc, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1ninc)

# Figure A9: Responsiveness spit by cantonal treatment and across electoral safety (incumbents 
# excluded) (export: 9 in X 4in)
# ------------------------------------------------------------------------------------------- #

m.beta <- fixef(M1ninc)
v.beta <- vcov(M1ninc)

# left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data_ninc$safety), 150, by = 1), 
             0*seq(min(data_ninc$safety), 150, by = 1),
             mean(data_ninc$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data_ninc$safety), 150, by = 1), 
             1*seq(min(data_ninc$safety), 150, by = 1),
             mean(data_ninc$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

plot (seq(min(data$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data$safety),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data$safety), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data$safety), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data$safety[data$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data$safety),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data$safety), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data$safety), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA9.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, 
   p.hat0, p.hat1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M1ninc)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A11 Robustness Check: Linear Probability Model

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A12: RI linear regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

M1lpm <- lmer(answer_r ~ 
               + treatment_party_r
               + treatment_cant_r
               + safety
               + I(treatment_cant_r*safety)
               + age 
               + male_r
               + party1 + party2 + party4 + party5 
               + lang_r
               + (1 | canton), data=data)
summary(M1lpm)

# Figure A10: Responsiveness spit by cantonal treatment and across electoral safety (linear 
# probability model) (export: 9 in X 4in)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

m.beta <- fixef(M1lpm)
v.beta <- vcov(M1lpm)

# left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data$safety), 150, by = 1), 
             0*seq(min(data$safety), 150, by = 1),
             mean(data$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.ci.XX0 <- apply(y.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data$safety), 150, by = 1), 
             1*seq(min(data$safety), 150, by = 1),
             mean(data$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.ci.XX1 <- apply(y.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

# left panel

plot (seq(min(data$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data$safety),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data$safety), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data$safety), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data$safety[data$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((y.hat1)-(y.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data$safety),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data$safety), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data$safety), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA10.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M1lpm)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A12 Robustness Check: Volunteer Effects

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A13: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

# m1

M1sender <- glmer(answer ~ 
               + treatment_cant_r
               + sender_male_r
               + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1sender)

# M1

M2sender <- glmer(answer ~ 
               treatment_cant_r
             + C(factor(data$sender))
             + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M2sender)

# clean up

rm(M1sender, M2sender)


# Table A14: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

M1valais <- glmer(answer ~ 
              treatment_party_r
            + treatment_cant_r
            + safety
            + I(treatment_cant_r*safety)
            + age 
            + male_r
            + party1 + party2 + party4 + party5 
            + lang_r
            + (1 | canton), data=data_valais, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1valais)


# Figure A11: Responsiveness spit by cantonal treatment and across electoral safety (Valais excluded)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

m.beta <- fixef(M1valais)
v.beta <- vcov(M1valais)

# left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data_valais$safety), 150, by = 1), 
             0*seq(min(data_valais$safety), 150, by = 1),
             mean(data_valais$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data_valais$safety), 150, by = 1), 
             1*seq(min(data_valais$safety), 150, by = 1),
             mean(data_valais$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

#left panel

plot (seq(min(data_valais$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data_valais$safety),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data_valais$safety), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data_valais$safety), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data_valais$safety[data_valais$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data_valais$safety),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data_valais$safety),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data_valais$safety), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data_valais$safety), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data_valais$safety), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA11.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, 
   p.hat0, p.hat1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M1valais)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# A13 Robustness Check: Candidates with Staff Excluded from Analysis

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table A15: RI logistic regression (outcome: answer)
# ------------------------------------------------------------------------------------------- #

M1staff <- glmer(answer ~ 
                  treatment_party_r
                + treatment_cant_r
                + safety
                + I(treatment_cant_r*safety)
                + age 
                + male_r
                + party1 + party2 + party4 + party5 
                + lang_r
                + (1 | canton), data=data_staff, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1staff)

# Figure A12: Responsiveness spit by cantonal treatment and across electoral safety (incum- bents excluded)
# ------------------------------------------------------------------------------------------- #

par(mfrow=c(1,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

m.beta <- fixef(M1staff)
v.beta <- vcov(M1staff)

# left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1), 
             0*seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1),
             mean(data_staff$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1), 
             1*seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1),
             mean(data_staff$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety, na.rm = TRUE),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data_staff$safety, na.rm = TRUE),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data_staff$safety[data_staff$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data_staff$safety, na.rm = TRUE),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data_staff$safety, na.rm = TRUE),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data_staff$safety, na.rm = TRUE), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# save chart

dev.copy(pdf,"figureA12.pdf", width=9, height=4)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, 
   p.hat0, p.hat1, XX0, XX1, y.hat0, y.hat1, m.beta, v.beta, M1staff)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# END

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

